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Abstract 


The upwind leapfrog or Linear Bicharacteristic Scheme (LBS) has previously been extended 
to treat lossy dielectric and magnetic materials. This paper examines different methodologies for 
treatment of the electric loss term in the Linear Bicharacteristic Scheme for computational electro- 
magnetics. Several different treatments of the electric loss term using the LBS are explored and 
compared on one-dimensional model problems involving reflection from lossy dielectric mate- 
rials on both uniform and nonuniform grids. Results using these LBS implementations are also 
compared with the FDTD method for convenience. 



1 Introduction 


Numerical solutions of the Euler equations in Computational Fluid Dynamics (CFD) have il- 
lustrated the importance of treating a hyperbolic system of partial differential equations with the 
theory of characteristics and in an upwind manner (as opposed to symmetrically in space). These 
two features provide the motivation to use the Linear Bicharacteristic Scheme (LBS), or the up- 
wind leapfrog method, for the construction of many practical wave propagation algorithms. The 
upwind leapfrog method has a more compact stencil compared with a classical leapfrog method. 
Clustering the stencil around the characteristic enables high accuracy to be achieved with a low 
operation count in a fully discrete way. This leads to a more natural treatment of outer bound- 
aries and material boundaries. The LBS treats the outer boundary condition naturally without 
nonreflecting approximations or matched layers. The interior point algorithm predicts the out- 
going characteristic variables, and the algorithm only requires information about the incoming 
characteristic variables at the domain boundaries. Through knowledge of the wave propagation 
angle, the local coordinates can be rotated to align with the characteristics, at which the boundary 
condition becomes almost exact. Therefore, no extraneous boundary condition or matched layers 
are required, which can introduce errors into the solution. The LBS also offers a natural treatment 
of dielectric interfaces, without any extrapolation or interpolation of fields or material properties 
near material discontinuities. 

The LBS was originally developed to improve unsteady solutions in computational acoustics 
and aeroacoustics [l]-[7]. It is a classical leapfrog algorithm, but is combined with upwind bias in 
the spatial derivatives. This approach preserves the time-reversibility of the leapfrog algorithm, 
which results in no dissipation, and it permits more flexibility by the ability to adopt a charac- 
teristic based method. The use of characteristic variables allows the LBS to treat the outer com- 
putational boundaries naturally using the exact compatibility equations. The LBS offers a central 
storage approach with lower dispersion than the Yee algorithm, plus it generalizes much easier to 
nonuniform grids. It has previously been applied to two and three-dimensional free-space elec- 
tromagnetic propagation and scattering problems [3], [6], [7]. It has also recently been extended 
to treat lossy dielectric and magnetic materials [8]. 

The objective of this paper is to examine different methodologies for treatment of the electric 
loss term in an attempt to find an accurate and self-consistent implementation that also works for 
perfect conductors in the limit of high conductivity (i.e. a — ► oo). The final goal is to develop 
an implementation that will allow for an accurate and efficient extension of this approach to two- 
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and three-dimensional problems. This is accomplished by examining several different implemen- 
tations of the electric loss term and by testing these implementations on one-dimensional model 
problems on both uniform and non-uniform grids. Section 2 of this paper provides a basic intro- 
duction to the LBS and it develops the different implementations for the electric loss term using 
the LBS. Section 3 presents results for one-dimensional model problems and Section 4 provides 
concluding remarks. 


2 Theory 


Maxwell's equations for linear, homogeneous and lossy media in the one-dimensional TE case 
(taking d/dy = d/dz = 0) are 
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dt 
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where a, a* are electric and magnetic conductivities, respectively. Using the electric displacement 
D = eE and making the substitution c = l/ v //7f gives 
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The idea for the LBS is to transform the dependent variables D y and H z to characteristic variables. 
The algorithm presented here is the simplest leapfrog scheme described by Iserles [9] combined 
with upwind bias, or simply, the Linear Bicharacteristic Scheme (LBS). To transform (3) and (4) 
into characteristic form, we first multiply (4) by c and then add and subtract from (3) to give 


d(D v + i c H t ) 

dt 

dt 


d(D y + \H z ) 
+ C dx 



+ —D y + a*ecH z = 0 
e y 

+ —D v — cr*ecH z - 0 
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Now define 


P = 



(7) 


to represent the left propagating solution and 

Q = D y + -H z 
c 


( 8 ) 
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to represent the right propagating solution. P and Q are otherwise known as the characteristic 
variables. Using these definitions, equations (5) and (6) can be rewritten as 


dQ dQ 1 fa 


W + ^ + 2{l~° CC 


(^_„. £C *)p + I(2 +17 v 2 )<3 = 0 


( 9 ) 

+ + ^)p + i(^-^)o = ° (10) 

It is convenient to define and store the following coefficients before the time-stepping begins 


(7 * 9 

a = — I - a ec 

e 

b = a ec 
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Equations (9) and (10) can be rewritten more simply as 

w + c 7r + l p + l e 3 = 0 

at ox 2 2 

OP dP a _ b„ n 

aT- c & + 2 P+ 2 Q = 0 


(11) 
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To develop the discretized algorithm for a one-dimensional system, the stencils of Figure 1 are 
proposed for the LBS. To solve the wave propagation problem without introducing dissipation, it 
is necessary that the stencil have central symmetry so that the scheme employed is reversible in 
time [2]. The stencil in Figure la is used for a right propagating wave and the stencil in Figure lb 
is used for a left propagating wave. The upwind bias nature of the stencils is thus clearly evident. 
References [1], [2], [5], [6] and [7] clearly show that the UL method is second-order accurate. Note 
that the last two terms in equations (13) and (14) represent the electric and magnetic loss (or source) 
terms. The goal of this paper is to determine the most accurate and efficient method of treating 
these terms. To that end, several different implementations for these source terms will now be 
developed. 


2.1 Method A 


The first method, denoted by Method A, is to index both loss terms at time level n + 1, for a 
semi-implicit formulation. This method was the original method recently proposed in [8]. Using 
the stencils shown in Figure 1, the resulting finite difference equations for (13) and (14) are 

+ C ( 0 " t jy+- + ^Q" + ' = 0 (15) 

- c ( ) + \Pi*' + jfer' - 0 (16) 
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Rearranging (15) and (16) gives 

Qi +l -Q? + Qi- 1 - Qtx + 2t/ (Q? - Q?_i) + 6AtTf +1 + aAtQ" +1 = 0 (17) 

p.n+1 _ pn + p, n+i _ p.n-1 _ 2|/ _ pnj + aAt pn+ 1 + & A fQ? +1 = 0 (18) 


where 1 / = cAt/Ax is the Courant number. Now, collecting common terms in (17) and putting the 
n and n — 1 time indexed terms on the right hand side (with a similar development for (18)) gives 

or 1 + (1 ^‘ A() ? +1 = [or, 1 + (1- 2v) (<3f -<?"_,)] /(1+aAt) (19) 

= [^+i i -(i-2-)w + i-^”)]/(i+«^) po) 

The right side of these equations is the residual and can be rewritten as 

Rx = [Qr-Ti 1 + (1 - 2i/) (Q? - Q?_,)] /(I + aAt) (21) 

SS = [P? + \ 1 -(l-2v)(P? +1 -Pn]/(l + aAt) (22) 

Now, in relation to equations (19) and (20), the following definitions are made 
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The update equations (19) and (20) are now rewritten as 


AX = R 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


(30) 


with a solution given by 


X = A~ X R 


where 


! 1 022 — O 12 
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(31) 
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The final update equations are then 


<?r +1 

= ~j («22 R\ ~ <*12 R 2 ) 
a 

(34) 

pn-f 1 
r i 

= ~\{-o. 2 \R 1 l ■}- a\\R 2 ) 
a 

(35) 


Note that no computational penalty is incurred by this matrix solution, since the a, b and matrix 
coefficients can be precomputed and stored before time stepping begins, thus increasing the effi- 
ciency of the method. Note that for a PEC as <r —* oo, from (34) and (35), Q" +1 = P™ +1 = 0, as 
required. For lossless (a = a* = 0), homogeneous dielectric and magnetic materials, the update 
equations are simply 


Q " +1 = (36) 

P ” +1 = Pi+i 1 ~ (1 — 2&0 (Pi+i — P?) (37) 


To implement the dielectric material interface boundary condition, consider the one-dimensional 
grid shown in Figure 2. The dielectric interface is located at grid point i, and the dielectric mate- 
rials can be lossy. The characteristic variables at grid point i, Pi and Qu are split into four com- 
ponents P\,Q\, P 2 and Q 2 . The terms Pi and Q\ exist just to the left of the material interface as 
shown in Figure 2. The remaining terms P 2 and Q 2 exist just to the right of the material interface. 
For material 1, equation (34) is used to predict the value for Q 7 ^ 1 at the boundary and for material 
2, equation (35) is used to predict the value for Pf + 1 at the boundary. Then the electromagnetic 
boundary conditions are used to solve for P^ 1 and Q 2 +1 in terms of the "known" characteristic 
variables Q " +1 and P 2 n+1 . To develop this procedure, the electromagnetic boundary conditions on 
the tangential field components are given by 

Ey\ = Ey 2 (3g) 

6l t 2 

H z 1 = H z2 (39) 

For the right-going wave, substituting (38) and (39) into (8) gives 

QT l = Dyt 1 + -KT 1 
^1 

= ^ (^ n+1 + (?2 +1 ) + ^ - P 2 +l ) (40) 

Similarly, substituting (38) and (39) into (7) yields 
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Since Q " +1 and P ^ +1 are known at the boundary point i, it is necessary to express Q^ 1 and 
in terms of these variables. Rearranging (40) and (41) gives 


where 


Qr 1 

pn+l 


= T 2 Q” +l 
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+ r 2 p 2 n+1 

+ r!P 2 n+1 
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Ti 
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_ 2 c 2 £j_\ 
c 2 e 2 + c ici/ 
c 2 e 2 - cieA 

c 2 e 2 + ciei J 

2 cie 2 \ 
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(42) 

(43) 

(44) 

(45) 

(46) 

(47) 


From (42) it is clear that the right-going wave in material 2 is a sum of a transmitted portion of 
the right-going wave in material 1 plus a reflected portion of the left-going wave in material 2. A 
similar argument can be made for the left-going wave in region 1. In fact, the coefficients in (42) 
and (43) are very similar to the classical Fresnel reflection and transmission coefficients. Because 
this is a semi-implicit method, a matrix solution is required to calculate P and Q at each grid point. 
Special care needs to be taken when the LBS calculates the solution at grid points i — 1 and i + 1 . 
At grid point i — 1 , the term P”_j in (22) becomes P". At grid point i, the terms Qf and P- 1 in (21) 
and (22) become Q" and P 2 n , respectively. At grid point i + 1 , the term Q”_j in (21) becomes Q 2 . 
Rearranging equation (17) we have for grid point i, 


(1 + aj At) Q7+ 1 = Q?Zi + (1 - 2u) (Q? - QU) - &iA tP? +1 


(48) 


Since the term P " +1 is unknown, we use the boundary condition (43) to give 
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A similar development using (18) and (42) yields 
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and we define 


ryn 

/l2 — 
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022 = 
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The matrix equation now becomes 
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The solution of this matrix equation is 

Q? +1 = ^(022^-012^) 
P? +1 = ^(-021^ + 011^) 


(54) 

(55) 

(56) 


(57) 


(58) 

(59) 


where d is defined by (33). Once Q" +1 and P% +1 have been obtained, the boundary conditions 
in (42) and (43) are applied to update Q% +1 and P" +1 . For the perfect conductor, the appropriate 
boundary conditions are Q% +1 = P 2 n+1 = 0.0 and P" +1 = — Q” +1 , where Q" +1 is predicted from 
(48). 


2.2 Method 6 


Method B averages the loss terms at time levels n and n+ 1, which still results in a semi-implicit 
formulation. The update equations remain the same as shown in (34) and (35), but the following 
definitions are made 


R’l 
. R* 


Qili + (1 — 2i/) (Ql-QL 1 ) 


aAt 


QP 


^]/( 1+ 


aAt\ 


2 

= ppj - (i - 2 u) (pt +1 - pp) - - ^pp] / (1 + 


(60) 

( 61 ) 


an = 
ai2 = 

1 

b At 

(62) 

(63) 

(2 “f* aAt) 
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(64) 

0-21 = 

(2 -j - aAt) 

0-22 = 

1 

(65) 


and the solution proceeds as in Method A. A dielectric surface boundary condition can be devel- 
oped for Method B in the same manner as that for Method A, but the details of that analysis are 
omitted for the sake of brevity. 
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2.3 Method C 


Method C indexes the Q loss term at time level n + 1 and the P loss term at time level n in 
(13). Similarly, in (14), the P loss term is indexed at n + 1 and the Q loss term is indexed at n. This 
method is an explicit approach and avoids the matrix solution inherent with the semi-implicit 
Methods A and B. To that end, the foilwing definitions are made: 


R? = 

Q7~i + (1 - 2i/) (Qi - Qti) - bAtPp 

(66) 

R% = 

P-+: \ l - (1 - 2*/) (PT +1 - PP) ~ bAtQ? 

(67) 

The update equations are then 

Q" +1 = R\/ (1 + aAt) 

(68) 


P^ +1 = P"/(l + aAt) 

(69) 


The dielectric surface boundary condition is obtained by direct application of (42) and (43). 


2.4 Method D 


Method D indexes both Q and P loss terms in (13) and (14) at time level n, which is also an 
explicit approach. The residual terms become 


Ri = Qi-I + (1 - 2 «/) (Qi - Qi-i) - aAtQft - bAPp 

R n 2 = Ptf -(! - 21/) (P? +1 - P?) - aAtPp - bAtQ* 
and the update equations are 


or 1 = Ri 



(70) 

(71) 


(72) 

(73) 


Unfortunately, this method is unstable for long time integrations and for high conductivity values. 
As noted by Thomas [5], this is a straightforward application of the upwind leapfrog method 
to the loss (or source) terms in (13) and (14), which he showed to be unstable. He proposed 
an alternative method for one source term which stabilizes the algorithm, and these methods 
including variations are considered next. Since Method D is unstable, no numerical results will be 
presented using this implementation. 


2.5 Method E 

Method E relies upon the method proposed by Thomas [5] to stabilize the upwind leapfrog 
algorithm when source terms are present in the governing equations. The idea is to introduce a 


9 


transformation 


(74) 


Q = e- at/2 V 


into (13) which gives (after some work) 


dV dV b ot/2 _ 

■W + C ^ + 2 e ' P 


= 0 


Applying a similar transformation 


p = e ~ atl2 w 


(75) 

(76) 


to (14) gives 


dt +c dx + 2 e Q 


(77) 


Method E and the following two methods deal with treatment of the remaining source term in 
(75) and (77). For Method E, we index the source terms at time level n + 1. Applying the stencils 
of Figure 1 to (75) and (77) and noting that 


yn _ e<inAt/2Qn 

W 71 — e anAt / 2 pp 


gives the following equations 

Q" +1 + bAtP? +1 = 

6AtQ" +1 + P" +1 = R% 

where 

R? = e~ aAt Q^Ii + e~ aAt/2 (l - 2v) (Qf - Q^_ x ) 
R% = -e- aAt ' 2 (l-2v)(P t n +1 -P % n ) 


(78) 

(79) 


(80) 

(81) 


(82) 

(83) 


Thomas [5] used a second order Taylor series expansion of the terms e _aAt and e aAt to derive 
the final update equations to avoid a costly exponential matrix evaluation and inversion. This is a 
concern especially for multidimensional applications and for multiple source terms. However, the 
Taylor series expansion is not implemented in the present context, since these exponential updat- 
ing coefficients can be precomputed and stored before time stepping begins. No computational 
penalty is incurred by using the full exponential terms. The update equations are 


Qr 1 = ^(«22 J*?-a 12 i$) 

ft 1 = i(-a 2 i m+anm) 


(84) 

(85) 
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with the following definitions 


an 

= 1 

(86) 

a 12 

= bAt 

(87) 

021 

= bAt 

(88) 

022 

= 1 

(89) 


and d is given by (33). The dielectric surface boundary condition is similar to that for Method A. 
Unfortunately, this approach exhibits a late-time instability problem for lossy dielectric materials. 
A Taylor series expansion of the exponential terms proposed by Thomas [5] results in a slightly 
different scheme that also exhibits late time instability. 


2.6 Method F 

Method F is similar to Method E, but averages the source term in (75) and (77) between time 
levels n + 1 and n. This results in the equations 


R n _ e -aAtQU-l +e -aAt/2^_ 

2u) (Q^Q tl )^^ e -aAt/2pn 

(90) 

R% = e~ aAt PJ!^ l - e" aA</2 (l - 

2v) (P? + i-P?)- b -^e- aAt / 2 Q? 

(91) 

Oil = 

1 

(92) 

012 = 

bAt/ 2 

(93) 

021 = 

bAt/2 

(94) 

022 = 

1 

(95) 


The final update equations are the same as (84) and (85). Again, the dielectric surface bound- 
ary condition is similar to that in Method A. This method also suffers from a late time instability 
problem, although the instability growth is not quite as rapid as in Method E. An alternative 
scheme based on a Taylor series expansion of the exponential terms also exhibits late-time insta- 
bility. Therefore, no results for Methods E and F will be presented. 


2.7 Method G 

The last method considered indexes the source term in (75) and (77) at time level n. This gives 
the final update equations as 

Q " +1 = e - a A«Qn-l +e -oAt/2(i_ 2l/ )(Qn_gn_ i ) _ 6Ate -aAi/ 2 p,n 

pn+\ = e -°*tpn-l_ e -aAt/2( 1 _ 2 „^pn +i _pn}_ hAte -aAt/2Qn 


11 


(96) 

(97) 


The dielectric surface boundary condition is obtained by a direct application of (42) and (43). This 
method is stable and does not exhibit any signs of late-time instability. 

3 Results 

The numerical results for this paper concentrate on one-dimensional model problems involv- 
ing reflection from lossy dielectric materials on both uniform and nonuniform grids. The problem 
space size is 1000 cells with nonperiodic boundary conditions. For the uniform grid, a space step 
of 1 cm is used, the time step is 0.33 ps and the Courant number v — 1. For the boundary condi- 
tions, a Gaussian point source at i = 0 is used to specify <2(0) and P(1000) — 0. For many complex 
geometries, it is often desirable to implement nonuniform grids to reduce the computational effort 
and memory resources and to improve modeling accuracy. We define a nonuniform grid by using 
a mesh stretch ratio of M = Ax max /Aa: m i n which is periodic every 10 cells. Figure 3 shows an 
expanded view of a typical one-dimensional nonuniform grid with a mesh stretch ratio of 2 and 
a periodicity of 10 cells. The Courant number for a nonuniform grid is defined by cAt/Ax m i n , 
where A x m i n is the smallest cell size in the nonuniform grid. 

The first problem is a reflection and transmission problem for a lossy dielectric half space on 
a nonuniform grid. The dielectric half space extends from 5 < x < 10 with material parameters 
e r2 = 10, cr 2 = 0.2, <j* = 0 and \x r2 — 1. Figure 4 shows the reflection coefficient magnitude re- 
sults for the exact solution, the FDTD method and the various LBS implementations calculated at 
x — 4.0 m (cell i — 276). Note that the FDTD results vary widely over the entire frequency band 
due to the nonuniform grid structure. However, most of the LBS implementations have compa- 
rable accuracy to the FDTD method, and these results are quite smooth, even for the nonuniform 
grid. Figure 5 shows the percent error in reflection coefficient magnitude for the LBS implemen- 
tations. The FDTD method was not shown due to the wide variation in error. Note that Method 
C for the LBS implementations has the overall lowest error over the entire frequency range, and is 
therefore the method of choice for modeling lossy media with the LBS. 

The second problem was the same as the previous case, but with a uniform grid of 1 cm spac- 
ing and material parameters e r2 = 50, a 2 — 0.2, a* = 0 and fi r2 = L Figure 6 shows the reflection 
coefficient magnitude results for the exact solution, the FDTD method and the various LBS imple- 
mentations calculated at x = 4.0 m (cell i = 401). Note that the FDTD method exhibits slightly 
better agreement for lower frequencies, but then diverges sharply at higher frequencies. Figure 
7 shows the percent error in reflection coefficient magnitude for the LBS implementations, which 
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maintain a rather flat frequency dependence at higher frequencies. This provides increased accu- 
racy and fidelity for the higher frequency components in a given simulation. Although Method 
B exhibits slightly lower error at very low frequencies for this particular problem, the differences 
between methods is not significant enough to prefer Method B over any other method. Method C 
is the most suitable alternative for multidimensional applications. 

The third problem was the same as the previous case, but with material parameters e r 2 = 80, 
<j 2 = 0.0, a* = 0 and \i r <i — 1. Figure 8 shows the reflection coefficient magnitude results for the 
exact solution, the FDTD method and the various LBS implementations calculated at x = 5.0 m 
(cell i = 500). Note that the FDTD method diverges sharply at higher frequencies, but the LBS 
methods are exact, and each LBS method reduces to the standard implementation for lossless di- 
electric media. 

The final problem was also the same, but with a lossless magnetic media with material param- 
eters e r 2 = 1, <72 = 0.0, cr* = 0 and fi T 2 = 5000. Figure 9 shows the reflection coefficient magnitude 
results for the exact solution, the FDTD method and the various LBS implementations calculated 
at x = 5.0 m (cell i — 500). Note how the FDTD method diverges from the exact solution, but 
the LBS methods are again exact for lossless magnetic media. Both the FDTD method and LBS 
methods had difficulty in predicting correct reflection coefficient results for lossy magnetic me- 
dia, which could be the result of not updating the electric and magnetic conduction currents as 
separate solution variables. This alternative will be explored for further research. 

4 Conclusions 

This paper has examined seven implementations of the Linear Bicharacteristic scheme for com- 
putational electromagnetics to treat heterogeneous lossy dielectric and magnetic materials. Two of 
the six implementations exhibited late-time instabilities and were therefore eliminated from fur- 
ther consideration. The remaining four implementations were tested on simple one-dimensional 
model problems involving reflection from lossless and lossy dielectric and/or magnetic half- 
spaces. It was found that the method of choice involved indexing one source term at the present 
time level, and the coupled source term at the previous time level. This implementation provided 
two benefits. First, it provided similar or better accuracy than the other implementations. Second, 
it is the most suitable implementation for multi-dimensional problems to avoid a matrix solu- 
tion at each time step. Overall, the LBS has several distinct advantages over conventional FDTD 
algorithms. First, the LBS is a second-order accurate algorithm which is about 2-3 times as eco- 
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nomical. The LBS can also be made to have zero dispersion error in certain instances. Second, 
the LBS provides a more natural and flexible way to implement surface boundary conditions and 
outer boundary conditions by using characteristics and an upwind bias technique popular in fluid 
dynamics. Third, the LBS provides more accurate results on nonuniform grids. The upwind bi- 
asing provides a more flexible generalization to unstructured grids. The results of this particular 
study and a previous study [8] indicate that the LBS is a superior algorithm for treatment of di- 
electric materials, especially its performance on nonuniform grids. Based on these results, the LBS 
is a very promising alternative to a conventional FDTD algorithm for many applications. 
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5 Figure Captions 

Figure 1. One-dimensional upwind leapfrog computational stencils for right-going (a) and left- 
going (b) characteristics. 

Figure 2. One-dimensional computational grid for the LBS showing characteristic variables, a di- 
electric interface located at cell i, and corresponding field components and characteristic variables 
used for the surface boundary condition. 

Figure 3. Section of a one-dimensional non-uniform grid with a mesh stretch ratio of 2 and a base 
cell size of 1 cm. The grid variation is periodic every 10 cells. 

Figure 5. Reflection coefficient magnitude versus frequency for reflection from a lossy dielectric 
half-space using the FDTD method and the LBS implementations on a non-uniform grid. 

Figure 6. Percentage error in reflection coefficient magnitude versus frequency for reflection from 
a lossy dielectric half-space using the LBS implementations on a non-uniform grid. 

Figure 7. Reflection coefficient magnitude versus frequency for reflection from a lossy dielectric 
half-space using FDTD and the LBS on a uniform grid. 

Figure 8. Percentage error in reflection coefficient magnitude versus frequency for reflection from 
a lossy dielectric half-space using the LBS implementations on a uniform grid. 

Figure 9. Reflection coefficient magnitude versus frequency for reflection from a lossless magnetic 
half-space using FDTD and the LBS on a uniform grid. 
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